1 Objectifs

L’objectif de ce travail est d’analyser la distribution statistiques et géographiques des teneurs en Carbone Organique (CO) de la BDAT pour les 5 périodes de temps suivantes : 1990-1994, 1995-1999, 2000-2004, 2005-2009 et 2010-2014. Dans un premier temps, le travail est porté à l’échelle de la France afin d’appréhender globalement les différences entre les périodes de temps. Le même travail est ensuite reconduit sur différentes strates géographiques pour affiner les résultats sur différentes échelles géographiques. Enfin, la troisième étape propose d’analyser la répartition spatiale des teneurs avec des facteurs explicatifs potentiels.

Les résultats présentés font suite à différents scripts de préparation de données dont la chaîne de traitements générale est consultable à cette adresse.

2 Analyse des teneurs en carbone organique par période (France)

Cette première étape a pour but d’analyser les teneurs en CO de la BDAT pour chacune des périodes de temps analysées. La cartographie des teneurs en CO est présentée dans un premier temps et les statistiques descriptives (courbes de fréquences cumulées et boxplot) sont présentées dans un second temps.

2.1 Représentation cartographique

2.1.1 Cartographie des effectifs de la BDAT

La cartographie du nombre d’analyse de la BDAT par canton pour les teneurs en CO est présentée dans la figure ci-dessous. Cette figure met en évidence une importante variabilité spatio-temporelle du nombre d’analyse. De façon générale (toutes périodes confondues), le nombre d’effectif est plus important sur la partie nord et ouest du pays. Le Massif Central et le Sud-Est sont faiblement concentrés en effectifs.

Les variabilité sont également marquées dans le temps. La période 9094 ne comporte pas d’analyse dans plusieurs régions comme la Normandie. Les autres périodes présentent une distribution des analyses plus stables.

Compte-tenu du manque de données observée dans la première période 9094, celle-ci ne sera pas prise en compte dans les analyses des teneurs en CO à l’échelle de la France.

2.1.2 Cartographie des teneurs en carbone organique (CO)

La cartographie des teneurs en CO de la BDAT est présentée dans les figures ci-dessous. Celle-ci est réalisée uniquement sur les 4 dernières périodes (9599, 0004, 0509, 1014) et appliquée sur un masque pour analyser les teneurs sur un nombre de canton commun entre les périodes.

La cartographie montre une distribution spatiale organisée et globalement similaire pour les différentes périodes analysées. De façon générale, cette organisation suit la lithologie du pays avec de fortes teneurs en carbone organique présentes dans les zones de socles et de piemonds et des valeurs plus faibles dans les principaux bassins sédimentaires (parisien et aquitain).

2.2 Statistiques descriptives

La figure 1 présente les courbes de fréquences cumulées des teneurs en carbone organique distribuées pour les 4 périodes identifiées. Les courbes de fréquences des 4 périodes présentent la même forme en “S” et s’individualisent juste avant le plateau, présentant une différence affectant les sols riches en teneurs organiques (entre 17 et 45 g/kg). Sur cette zone (figure à droite), la figure montre un décalage des courbes des périodes 2000-2004, 2005-2009 et 2010-2014 vers des valeurs plus faibles. Parmi ces 3 périodes, la période 2005-2009 est celle qui se décale le plus vers des teneurs plus faible tandis que la période 2010-2014 se rapproche des valeurs de 1995-1999, présentant ainsi une inversion de la tendance observée.

Ces observations mettent en évidence une diminution des teneurs en carbone entre les périodes 1990-1999 et 2000-2009 et une légère augmentation pour la période 2010-2014.

La distribution des teneurs en carbone organique par période est présentée figure 2 et les statistiques associées dans le tableau ci-dessous. La tendance de diminution des teneurs observée dans la figure 1 est également constatée dans ces deux éléments. La période 2000-2004 montre la valeur médiane la plus faible avec une valeur de 13.52 g/kg. Les valeurs les plus importantes sont observées pour les périodes 1995-1999 et 2010-2014 avec des médianes de teneurs en carbone organique équivalente, à 14 g/kg. Ces évolutions sont très légèrement marquées sur la figure 2 où l’évolution moyenne des teneurs en carbone organique baisse légèrement après la période 1995-1999.

Statistiques descriptives des teneurs en carbone organique pour les différentes périodes
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	5.23	11.05	14	15.51	18	82.36
0004	4.9	10.67	13.52	15.05	17.44	90.23
0509	4.1	10.8	13.7	15.12	17.57	69.2
1014	5.2	11	14	15.33	17.8	53.2

Les résultats du test de Wilcoxon présentés ci-dessous montrent que les différences des médianes cantonales à l’échelle de la France entre les périodes sont significatives pour la baisse observée entre les périodes 1995-1999 et 2000-2004, 1995-1999 et 2005-2009 et l’augmentation des teneurs entre les périodes [2000-2004 et 2010-2014]. Ces résulats confirment les tendances observées sur les courbes de fréquences cumulées (voir figure 1).

pairwise.wilcox.test(melted.bdat[,"value"], melted.bdat[,"annees"])

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.bdat[, "value"] and melted.bdat[, "annees"] 
## 
##      9599   0004   0509  
## 0004 0.0017 -      -     
## 0509 0.0348 0.5975 -     
## 1014 0.7060 0.0043 0.0678
## 
## P value adjustment method: holm

2.3 Conclusion

La médiane cantonale des teneurs en CO des analyses de la BDAT a évolué significativement entre les périodes étudiées (de 1995 à 2014). De 1995 à 2005-2009, la médiane des médianes des teneurs a diminué, passant de 14 g/kg à 13.70 g/kg. Cette tendance s’est inversée significativement entre les périodes 2000-2000 et 2010-2014 pour finalement revenir aux teneurs initiales de 1995-1999, à 14 g/kg. Au plus fort de ces évolutions, la baisse médiane des médianes des teneurs en CO entre 1995-1999 et 2000-2004 est de 3.4%, soit environ une baisse annuelle de 0.7% par an. A l’inverse, on retrouve un taux d’augmentation plus faible entre les périodes 2000-2004 et 2010-2014 avec 0.3% par an.

Ces résultats sont comparables avec les travaux de Bellamy et al., 2005. Ces auteurs ont mis en évidence sur un ensemble d’analyses de sol réalisé en Angleterre et au Pays-de-Galles de 1978 à 2003 une baisse moyenne de teneurs en CO de 0.6% par an (relatif par rapport aux teneurs initiales).

Ce travail d’analyse sur l’ensemble du territoire est reconduit sur plusieurs strates géographiques afin d’étudier les évolutions sur des compartiments environnementaux homogènes, du point de vue du climat et de l’occupation du sol.

3 Analyse des teneurs en carbone organique par période et par différentes strates géographiques

L’analyse des teneurs en CO est appliquée à l’échelle de plusieurs de strates. Ces strates représentent des entités géographiques homogènes d’un point de vue climatique, du type d’élevage. L’analyse par strate peut permettre d’observer des tendances régionales à l’intérieur desquelles les mécanismes d’évolution des teneurs en CO sont supposés homogènes.

3.1 Analyse par type de climat

Compte tenu de l’importance des variables climatiques dans le contrôle de la distribution spatiale des teneurs en carbone organique, une stratification des analyses précédemment réalisées est proposée. Celle-ci est basée sur une typologie des types de climat Joly et al., 2010. La typologie est consultable sur cette page.

La construction des données climatiques est consultable sur ce fichier FS_bdd_elab_climat.Rmd

Le tableau ci-dessus présente le nombre de canton pourvu en analyse de la BDAT par type de climat. Les types de climats 6, 7 et 8 englobent moins de 100 cantons. Ils seront écartés de l’analyse pour des raisons Parmi les strates restantes, la variabilité dans le nombre d’analyse est importante, passant de 163 pour le type 1 à 722 pour le type 3. Ces différences seront à prendre en compte dans l’analyse des résultats.

Nombre d’analyse par type de climat
Type_climat	Nbr
1	163
2	274
3	722
4	440
5	374
6	32
7	78
8	95

3.1.1 Cartographie des types de climat

La carte ci-dessous présente la typologie proposée. Pour une description des types de climat, voir Joly et al., 2010.

3.1.2 Statistiques descriptives

On commence par analyser les teneurs par strates géographiques… Ci-dessous, la figure présente la répartition des teneurs en C pour les différentes périodes et pour les principales régions d’élevage. Rajouter les différences significatives

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.bdat_clim[, "value"] and melted.bdat_clim[, "typo_clim"] 
## 
##   1       2       3       4      
## 2 < 2e-16 -       -       -      
## 3 < 2e-16 < 2e-16 -       -      
## 4 < 2e-16 < 2e-16 < 2e-16 -      
## 5 < 2e-16 2.6e-10 < 2e-16 < 2e-16
## 
## P value adjustment method: holm

La figure 4 présente les courbes de fréquence cumulées pour les différents types de climats. Les courbes de fréquences cumulées pour les 5 types de climats sélectionnés présentent une forme similaire. Les courbes des différentes périodes tendent vers des valeurs de teneurs en CO plus faible sur une large gamme de teneur pour les types de climat 1, 2, 3 et 4. Dans ces 2 cas, les périodes 0004 et 0509 présentent des valeurs plus faibles que les teneurs 9599 et 1014. Dans le type de climat 4, l’évolution des teneurs s’accentue sur la gamme 17-30 g/kg. Sur cette gamme, la période 1014 est particulièrement contributrice de la baisse des teneurs. A l’inverse, la période 0509 montrent des teneurs en CO plus forte pour cette gamme de valeur. Dans la zone climatique 2, les périodes 0004 et 0509 présentent des teneurs en CO plus faible que les teneurs initiales et celles de la période 2010-2014. Ces évolutions impliquent l’ensemble de la gamme des teneurs en CO.

La courbe de fréquence cumulée pour le type de climat 5 présente une dynamique similaire avec des changements important pour les teneurs en CO de 20-30 g/kg. Dans cette situation, la période 1014 présente également une courbe avec des teneurs en CO plus faible que la période 9599 et des teneurs en CO plus forte pour 0509 et 0004.

Cette figure montre des tendances comparables entre les types de climats. L’évolution des teneurs en CO est différente en fonction des gammes des teneurs en CO : la période 1014 présente une forte diminution sur les teneurs 17-30 g/kg et une augmentation sur les teneurs inférieures. Les périodes 0004 et 0509 affichent une tendance inverse.

Statistiques descriptives des teneurs en carbone organique pour le type de climat 1. Nombre de canton analysé : 163
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	10	17	20.58	24.13	27	82.36
0004	10.85	16.96	20.06	24.03	25.4	90.23
0509	8.7	16.59	20.2	23.25	24.12	69.2
1014	12	17.02	20.76	23.77	26.5	53.2

Statistiques descriptives des teneurs en carbone organique pour le type de climat 2. Nombre de canton analysé : 274
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	9	14	16.96	17.35	20	34.8
0004	7.85	13	15.68	16.84	19	41.18
0509	8.6	13.12	15.75	16.62	18.6	44.74
1014	10.5	14	16.3	17.2	19.2	37

Statistiques descriptives des teneurs en carbone organique pour le type de climat 3. Nombre de canton analysé : 722
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	5.35	10.54	12.27	13.14	15.41	30.35
0004	4.9	10.12	11.65	12.66	14.77	33.37
0509	5.12	10.44	11.7	12.73	14.66	28.14
1014	6.38	10.88	12.2	13.09	14.82	32.8

Statistiques descriptives des teneurs en carbone organique pour le type de climat 4. Nombre de canton analysé : 440
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	5.81	11.6	14.04	14.53	16.71	37
0004	5	11	13.78	13.9	16.18	36.86
0509	5.84	10.88	13.94	14.16	16.9	30.2
1014	6.29	11.29	14.18	14.33	16.71	34.8

Statistiques descriptives des teneurs en carbone organique pour le type de climat 5. Nombre de canton analysé : 374
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
9599	8.14	13.61	17.96	19.12	23.65	40.3
0004	7.77	13.08	17.84	18.71	23.12	40.12
0509	6.8	13.74	18	19.04	23.22	45
1014	8.39	13.56	17.75	18.71	23.15	42

La figure 5 représente la distribution des teneurs en carbone organique sous forme de boxplot pour les différentes périodes et les différentes zones climatiques. Les tableaux ci-dessus présentent les statistiques descriptives des teneurs en CO pour les différentes périodes et les différents types de climat. De ces tableaux et de ces boxplots, les zones 1, 2, 3 et 4 présentent une dynamique similaire avec une baisse puis une augmentation des teneurs. Ces évolutions se rapprochent des tendances nationales précédemment observées. Entre ces zones, seul les résultats pour le type de climat 3 sont significatifs. La baisse des teneurs entre 9599 et 0004 et l’augmentation des teneurs en CO entre 0004 et 1014 et 0509-1014 sont significatives.

La zone 5 présente une évolution en dent de scie, avec une augmentation, une diminution et une augmentation. Ces évolutions ne sont pas significatives d’après le test de Wilcoxon.

## $`1`
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.test[, "value"] and melted.test[, "annees"] 
## 
##      9599 0004 0509
## 0004 1    -    -   
## 0509 1    1    -   
## 1014 1    1    1   
## 
## P value adjustment method: holm 
## 
## $`2`
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.test[, "value"] and melted.test[, "annees"] 
## 
##      9599 0004 0509
## 0004 0.2  -    -   
## 0509 0.1  1.0  -   
## 1014 1.0  0.3  0.2 
## 
## P value adjustment method: holm 
## 
## $`3`
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.test[, "value"] and melted.test[, "annees"] 
## 
##      9599    0004    0509   
## 0004 0.01661 -       -      
## 0509 0.08612 0.68338 -      
## 1014 0.68338 0.00084 0.01661
## 
## P value adjustment method: holm 
## 
## $`4`
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.test[, "value"] and melted.test[, "annees"] 
## 
##      9599 0004 0509
## 0004 0.19 -    -   
## 0509 0.89 1.00 -   
## 1014 1.00 0.30 1.00
## 
## P value adjustment method: holm 
## 
## $`5`
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.test[, "value"] and melted.test[, "annees"] 
## 
##      9599 0004 0509
## 0004 1    -    -   
## 0509 1    1    -   
## 1014 1    1    1   
## 
## P value adjustment method: holm

3.2 Analyse par région d’élevage

Les grandes régions d’élevage représentent un zonage des principaux systèmes d’élevage en France. Le zonage a été construit sur la base de critères pédo-climatiques et d’un regroupement des petites régions agricoles. Par construction, ce zonage propose donc un regroupement de zones où les conditions d’apport, de stockage et d’évolution de carbone dans le sol peuvent être homogène.

Il représente donc un zonage intéressant à exploiter pour étudier les évolutions des teneurs en CO à l’échelle de la France sous le regard de compartiments fonctionnels homogènes.

Nombre de canton possédant des analyses par type de région d’élevage
Region_elevage	Nbr
A	616
B1	436
B2	168
C1	274
C2	76
D	163
E1	139
E2	136
F1	11
F2	64
G1	42
G2	37
H	16

Le tableau ci-dessus présente le nombre de canton ayant des analyses en teneurs en CO par type de région d’élevage. Les régions d’élevages où le nombre de canton est inférieur à 100 sont éliminées. Il reste donc les régions suivants :

A : Zones de grandes cultures ou sans élevage
B1 : Zone de polyculture-élevage du Bassin Aquitain & Rhône-Alpes & Alsace
B2 : Zone de polyculture-élevage du Bassin Parisien
C1 : Zone intensive du Grand Ouest (zone laitière avec alternatives à l’élevage)
D : Zone herbagère du Nord-Ouest
E1 : Zone herbagère du Nord-Est (de tradition laitière)
E2 : Zone herbagère du Nord Massif-Central (de tradition allaitante)

Parmi les strates restantes, la variabilité dans le nombre d’analyse est importante, passant de 163 pour le type 1 à 722 pour le type 3.

3.2.1 Cartographie

3.2.2 Statistiques descriptives

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  melted.bdat_regelevage[, "value"] and melted.bdat_regelevage[, "zonage_cplt"] 
## 
##    A       B1      B2      C1      D       E1   
## B1 < 2e-16 -       -       -       -       -    
## B2 0.796   < 2e-16 -       -       -       -    
## C1 < 2e-16 < 2e-16 < 2e-16 -       -       -    
## D  < 2e-16 < 2e-16 < 2e-16 0.796   -       -    
## E1 < 2e-16 < 2e-16 < 2e-16 1.3e-05 3.1e-05 -    
## E2 < 2e-16 < 2e-16 < 2e-16 5.9e-12 1.4e-11 0.015
## 
## P value adjustment method: holm

La figure 7 présente les courbes de fréquences cumulées pour les différentes régions d’élevage. Sur cette figure, on remarque l’effet les différences de la densité d’analyse par région. Les régions C2, G1 et G2 présentent des courbes en marche d’escalier, illustrant un nombre insuffisant de données. Mise à part ces différences, la forme des courbes est globalement similaire pour les différentes strates des régions d’élevage. En revanche, les tendances d’évolution des teneurs en CO sont différentes. Par exemple, dans la région C2, la fréquences des teneurs fortes en CO augmente pour les périodes 0004-0509 tandis que ces mêmes périodes affichent des diminutions importantes sur les zones D et E1 notamment.

A ce niveau, il est difficile de rentrer dans le détail. D’autres figures sont produites ci-dessous pour améliorer l’analyse de l’évolution des teneurs au niveau des différentes régions d’élevage.

Avant tout, les tables ci-dessous présentent les résultats des tests de Wilcoxon. Mis à part les zones d’élevage A et dans une moindre mesure B1, l’ensemble des régions d’élevage ne présente pas de différences significatives. Attention, il s’agit d’analyse à l’échelle de la région d’élevage et non du canton.

Ces résultats atténuent les observations visibles sur les figures ci-dessous (voir figure 8). Celles-ci présentent des évolutions frappantes et différentes selon les régions, mais elles ne sont pas significatives.

3.3 Conclusion

L’étude des teneurs en CO au niveau des différentes strates géographiques de type climatique et d’élevage a revélé des difficultés d’analyse en raison d’un nombre limité de données par strates. Les tendances observées à l’échelle des strates ne sont pas significatives et le faible nombre d’analyse dans les strates peut expliquer ces résultats.

Toutefois, sur l’ensemble des périodes analysées, les différences des teneurs entre les différentes strates géographiques sont systématiquement différentes. Ce qui démonter l’intérêt de travailler sur ce niveau d’échelle.?

4 Analyse des facteurs contrôlant la distribution spatiale

Dans cette partie, les facteurs contrôlant la distribution spatiale des teneurs en carbone organique des différentes périodes sont analysés. Le travail est dans un premier temps porté par une analyse en composante principale pour identifier les relations entre les facteurs potentiellement explicatifs et identifier les principales variables d’intérêts. Dans un deuxième temps par une modélisation avec la méthode des arbres de régression boostés pour affiner l’analyser et classer les variables explicatives par ordre d’importance (contribution).

4.1 Analyse en composante principale

La figure 9 présente la distribution des variables sur les axes 1 et 2 de l’ACP. Environ 60 pourcent de l’information est contenu dans ces deux premiers axes. Description de la distribution des variables dans l’espace des ACP :

L’axe 1 représente 42% de l’information. Cet axe est fortement associé à l’occupation du sol. Les différentes variables associées aux grandes cultures et à l’élevage sont bien représentées et ont une forte contribution dans la variance du jeu de données. Les statistiques liées à la part d’STH, de prairies, de SFP et d’OTEX polyculture élevage sont corrélées et sont opposées aux variables associées aux grandes cultures (proportion des surfaces de céréales/SAU et OTEX grandes cultures).
L’axe 2 représente 18% de l’information. La correlation avec les variables est moins nette que l’axe 1. Les variables de type climat et topographie sont légèrement correlées à cet axe. D’autres variables d’occupation du sol sont associées à l’axe 2. Celles-ci concernent les zones agricoles hétérogènes (\(clc_25_90\)), les zones de maïs fourrage et les zones qui ont une densité importante en UGBTA. Logiquement, ces deux dernière variables sont liées.
Les variables liées à Corine Land Cover ont un faible impact dans l’espace de corrélation. Etant donné qu’elle représente une information similaire aux données du recencement agricole, elle sont éliminées du jeu d’analyse par la suite.*

Dans cet espace, la teneur en carbone organique initiale (période 90-94) est moyennement représentée (-50% de contribution). Cette variable est directement opposée aux variables climatiques ce qui met en valeur l’importance de ce facteur dans la distribution spatiale des teneurs à l’échelle de la France. Dans une moindre mesure, la teneur en CO est correlée aux surfaces fourragères principales et aux UGBTA.

4.2 Test d’endogénéité des variables d’occupation du sol

L’objectif de cette partie est d’évaluer l’endogénéité des différentes variables d’occupation du sol. La notion d’endogénéité provient de l’économétrie et résulte de la corrélation entre un paramètre étudié et son résidu. Dans notre situation, l’hypothèse d’endogénéité est suspectée sur les variables d’occupation du sol car celle-ci influencent les teneurs en carbonique mais peuvent également dépendre de la teneur en CO du sol.

Le code ci-dessous présente l’application du test de Wu-Hausman. L’application de ce test demande l’utilisation d’instrument pour tester l’endogénéité. Les teneurs en CO pour des périodes postérieures aux périodes de l’occupation du sol ont été utilisées.

Le tableau ci-dessous présente les résultats du test appliqué pour toutes les variables d’occupation du sol. Le test de Sargan permet de tester la validité des instruments utilisés et les p_value du test d’Hausman permettent d’identifier les variables endogènes.

Les p-value sont supérieures à 0.05 et l’hypothèse null ne peut être rejeté. L’ensemble des variables testées ne sont pas endogènes. Les intruments utilisés sont tous valides à l’exception du test appliquée sur l’OTEX grandes cultures en 1988.

# selon http://eclr.humanities.manchester.ac.uk/index.php/IV_in_R
library(AER)

## Loading required package: car

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

# Ensemble des variables d'occupation du sol suspectées d'être endogène
endo <- c("p_sfp1988","elevage1988","elevagehorsol1988","polyelevage1988","p_prairie1970","p_prairie1979","p_prairie1988","p_sfp1970","p_sfp1979","p_mf1988","p_sth1970","p_sth1979","p_sth1988","p_mf1970","p_mf1979","p_mf1988","ugbta1988","grdcultures1988","p_cop1970","p_cop1979","p_cop1988","p_c1988","p_c1970","p_c1979")

# Ensemble des variables exogènes (climat+topo+pédo)
exo <- c("altimean","hpluie_an","jchauds_an","jfroids_an","pluie_ecart_janv","std_pluie_juil","std_temp_janv","std_temp_juil","argi_med","sabt_med")
exoplus <- paste(exo,collapse="+")
variamodel <- c(endo,exo)

# Instrument (teneur en carbone sur deux périodes postérieur aux occupations du sol(>2000))
instru <- c("corgox1014","corgox0509")
instruplus <- paste(instru,collapse="+")

vNames <- unique(c("corgox9599",instru,variamodel))
d <- dcast.bdat[complete.cases(dcast.bdat[vNames]),vNames]

# Test des instruments et de l'endogénéité
intrutest <- list()
HausWu <- list()
for(i in endo){
  varia <- variamodel[!variamodel %in% i]
  
  # Régression linéaire entre la variable endogène et les instrus et les variables exogènes
  first_stage <- lm(paste(i, " ~ ", paste(c(exoplus,instruplus),collapse="+"),sep=""),data=d)
  
  instrHausWutest <- waldtest(first_stage,.~.-corgox1014-corgox0509)
  intrutest[[i]] <- instrHausWutest[2,4]#Sargan

  # Régression linéaire du carbone (1995-1999) avec les variables exogène et le résidus de la première régression linéaire
  Hausman_reg <- lm(paste(vNames[1], " ~ ", paste(exo,collapse="+"),"+first_stage$residuals",sep=""),data=d)#Ajout des résidus au modèle de base
  
  # Comparaison des deux modèles 
  HausWutest <- waldtest(Hausman_reg,.~.-first_stage$residuals)
  HausWu[[i]] <- HausWutest[2,4]
}

rest <- cbind.data.frame(melt(intrutest),melt(HausWu))[,1:3]
colnames(rest) <- c("Sargan_test","Variables","Hausman_test")
pander(rest,caption="Résutlats des tests d'endogénéité (test de Wu-Hausman)")

Résutlats des tests d’endogénéité (test de Wu-Hausman)
Sargan_test	Variables	Hausman_test
1.939e-07	p_sfp1988	0.4638
1.671e-11	elevage1988	0.5918
0.02057	elevagehorsol1988	0.7793
4.728e-07	polyelevage1988	0.5593
2.546e-05	p_prairie1970	0.4355
9.883e-06	p_prairie1979	0.3204
3.184e-06	p_prairie1988	0.4167
1.48e-06	p_sfp1970	0.423
3.267e-07	p_sfp1979	0.421
4.098e-09	p_mf1988	0.8531
0.002191	p_sth1970	0.2216
0.0004009	p_sth1979	0.2028
0.0001678	p_sth1988	0.1869
2.624e-11	p_mf1970	0.8579
1.705e-11	p_mf1979	0.914
0.0002275	ugbta1988	0.4242
0.125	grdcultures1988	0.7444
0.006611	p_cop1970	0.3454
0.003452	p_cop1979	0.3642
0.01577	p_cop1988	0.4096
0.0003825	p_c1988	0.6991
0.003479	p_c1970	0.3521
0.001097	p_c1979	0.3105

4.3 Sélection des variables d’intérêt

L’objectif de cette partie est de sélectionner les principales variables explicatives et d’étudier l’effet des variables d’origines anthropiques sur un modèle linéaire des teneurs en CO. Pour cela, l’algorithmme suivant est appliqué :

1 : Pour chaque type de variables explicatives (naturelle et anthropique), application d’un modèle linéaire avec sélection du meilleur modèle sur le critère AIC. La méthode utilisée pour construire les modèles est basé sur le stepwise/backward (fonction stepAIC).
2 : Extraction des meilleurs modèles et calcul du R2 non-ajusté
3 : Construction d’un modèle linéaire basé sur les variables anthropiques et naturelles précédemment sélectionnés. Calcul du R2 non-ajusté et comparaison avec les modèles linéaires anthropiques et naturelles.

Les résultats de ces tests sont présentés dans le tableau ci-dessous. Pour l’ensemble des teneurs en CO modélisées, la part des variables de types anthropiques est toujours la plus faible. Elle se situe entre 5 et 7% pour les (voir deuxième tableaux ci-dessous). La part des variables naturelle (topographie, pédologie et climat) est bien plus importante.

Teneur_CO	Type_variable	R2_nonajusté
corgox9599	Anthropique	0.3628
corgox9599	Naturelle	0.5345
corgox9599	Complet	0.5808
corgox0004	Anthropique	0.3861
corgox0004	Naturelle	0.5522
corgox0004	Complet	0.6181
corgox0509	Anthropique	0.401
corgox0509	Naturelle	0.536
corgox0509	Complet	0.6052
corgox1014	Anthropique	0.3457
corgox1014	Naturelle	0.5287
corgox1014	Complet	0.5751

Contribution (en %) des variables anthropiques dans la modélisation des teneurs en carbone organique pour les différentes périodes
corgox9599	corgox0004	corgox0509	corgox1014
4.63	6.58	6.93	4.64

4.4 Régressions linéaires

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name              grob
## 1 1 (1-1,1-1) arrange   gtable[arrange]
## 2 2 (2-2,1-1) arrange gtable[guide-box]

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name              grob
## 1 1 (1-1,1-1) arrange   gtable[arrange]
## 2 2 (2-2,1-1) arrange gtable[guide-box]

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name              grob
## 1 1 (1-1,1-1) arrange   gtable[arrange]
## 2 2 (2-2,1-1) arrange gtable[guide-box]

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name              grob
## 1 1 (1-1,1-1) arrange   gtable[arrange]
## 2 2 (2-2,1-1) arrange gtable[guide-box]

4.5 Modèle de fouille de données

Dans cette partie, on testera rapidement deux modèles de fouille de données utilisés dans plusieurs travaux sur les sols (voir Murciano et al., 2015 pour une description des modèles)

GBM : Le principe de cet algorithme est de construire une séquence de modèles afin qu’à chaque étape, chaque modèle ajouté à la combinaison apparaisse comme un pas vers une meilleure solution.
Cubist : Cubist est un modèle de régression axé prédiction qui combine les idées de modèle à base de règles, décrites dans Quinlan (1992) avec des corrélations supplémentaires basées sur les voisins les plus proches dans l’ensemble d’apprentissage (voir Quinlan (1993) pour plus de détails).L’algorithme de Cubist est basé sur des règles permettant l’accumulation de modèles prédictifs basée sur des partitions du jeu de données.

4.5.1 Comparaison des modèles (GBM et Cubist)

Ci-dessous, la figure et les tableaux présentent les résultats de la qualité de la prédicition des deux modèles sur les 10 validations croisées. Le modèle cubist présente de meilleurs résultats que le modèle gbm. Pour la précision du modèle, la médiane de la racine quarré de l’erreur quadratique moyenne du modèle cubist est de 2.1384932 tandis qu’elle est de 2.205384 pour gbm. Les coefficients de détermination entre les deux modèles sont quasiment identique, avec un légèr avantage pour le modèle cubist. Le modèle Cubist sera utilisé par la suite.

## 
## 
##   * **RMSE**:
## 
##     ---------------------------------------------------------------------
##        &nbsp;     Min.   1st Qu.   Median   Mean   3rd Qu.   Max.   NA's
##     ------------ ------ --------- -------- ------ --------- ------ ------
##       **gbm**    2.033    2.169    2.205    2.29    2.456   2.699    0
## 
##      **cubist**   1.89    2.068    2.138   2.168    2.193   2.637    0
##     ---------------------------------------------------------------------
## 
##   * **Rsquared**:
## 
##     ---------------------------------------------------------------------
##        &nbsp;     Min.   1st Qu.   Median   Mean   3rd Qu.   Max.   NA's
##     ------------ ------ --------- -------- ------ --------- ------ ------
##       **gbm**    0.7358  0.8197    0.8302  0.8226   0.851   0.8722   0
## 
##      **cubist**  0.7392  0.8163    0.8344  0.8333  0.8722   0.891    0
##     ---------------------------------------------------------------------
## 
## 
## <!-- end of list -->

4.5.2 Application du meilleur modèle Cubist

Ci-dessous, la modélisation avec Cubist est relancé selon l’algorithme suivant :

1. Définition des paramètres de modélisation
   - Validation croisée sur les différents paramètres de modélisation (10 fois)
   - Sélection des meilleurs paramètres
2. Application du modèle avec les meilleurs paramètres (100 fois)
   - Sélection aléatoire de 80% du jeu de données pour l'apprentissage et 20% pour la validation
   - Construction du modèle sur 80% du jeu de données
   - Calcul de l'importance des variables dans le modèle
   - Prédiction sur 20% du jeu de données et comparaison avec le jeu de données initial
   - Calcul des indicateurs de qualité (R2, MSE, RMSE)
3. Calcul de la moyenne des indicateurs sur les 100 répétitions
   - Calcul de la moyenne des indicateurs de qualité sur le nombre de répétition
   - Calcul de l''importance moyenne des variables pour le nombre de répétition

Les résultats de la validation croisée de la modélisation des teneurs en CO pour les différentes périodes sont présentés dans le tableau ci-dessous. Sur les 100 répétitions, l’application du modèle sur les 20% du jeu d’apprentissage montre une capacité de prédiction des teneurs en CO équivalente entre les périodes et de bonne qualité. La modélisation des teneurs en CO pour la période 2010-2014 montrent la capacité de prédiction la plus mauvaise avec un coefficient de détermination (R2) de 0.78 et un erreur moyenne (Root Mean Square Error) de 2.4 g/kg. La période 2000-2004 affiche les meilleurs avec un R2 de 0.83 et une erreur moyenne de 2.2 g/kg. En moyenne la modélisation des teneurs en CO toute période confondues affiche donc une capacité d’explication de la variance de plus de 80%.

Moyenne des indicateurs de la validation croisée pour les différents modèles
	r2	MSE	RMSE
corgox9599	0.8235	5.199	2.277
corgox0004	0.8351	5.12	2.259
corgox0509	0.8173	5.701	2.385
corgox1014	0.7841	6.015	2.445

Les 15 variables explicatives de ces modélisations sont présentées dans la figure ci-dessous. Le taux d’argile (argi_med) et l’écart de pluie à la moyenne en janvier (pluie_ecart_janv) sont les deux premières variables explicatives des 4 modèles. Elles présentent toutes deux un pourcentage d’importance supérieur à 50%. L’altitude moyenne (altimean) par canton et le pourcentage de surface fourragère principale (pour différentes périodes) (p_sfp) sont également des variables explicatives d’importance.

Cette dernière variable est quasiment présente sous toutes ces déclinaisons temporelles (1970, 1979, 1988, 2000, 2010). Dans une moindre mesure, le pourcentage de prairie dans la SAU et le pourcentage de maïs fourrage sont les autres variables de type occupation du sol qui ressortent régulièrement dans ce classement des 15 variables explicatives les plus importantes.

L’ordre d’importance des variables explicatives des 4 périodes des teneurs en CO ne sont pas complètement identique à défaut de la première variables explicatives (l’argile). Toutefois, on note de nombreuses similitude entre les 3 premiers modèles (1995-1999 à 2005-2009). La dernière modélisation (2010-2014) présente des résultats différents, avec une part de variable d’origine climatique plus importante.

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name              grob
## 1 1 (1-1,1-1) arrange   gtable[arrange]
## 2 2 (2-2,1-1) arrange gtable[guide-box]

En conclusion Ces résultats montrent que la répartition des teneurs en CO de la BDAT à l’échelle de la France est principalement expliquée par des variable d’origine naturelle (pédologique, climatique, et topographique). L’effet de l’occupation du sol apparaît secondaire. Ces résultats concordent avec les premières conclusions du travail de la modélisation avec des régressions linéaires multiples effectuées avec des co-variables naturelles ou d’origines anthropiques.

4.6 Conclusion

A la suite de ce premier travail, la modélisation des teneurs en CO doit être reconduite à l’échelle d’un zonage où le contexte climatique, pédologique et topographique est homogène afin d’identifier les variables d’occupation du sol susceptibles d’agir sur les teneurs en CO.

5 Analyse par stratification

Cette section reprend le travail de modélisation réalisé précédemment en l’appliquant à l’échelle des strates liée au type de climat et aux principales régions d’élevage.

5.1 Analyse linéaire par stratification

Attention, dans les résultats, revoir le calcul car valeur R2 abérante. Prendre en compte la différences dans le nombre d’échantillon par strates.

5.2 Analyse Cubist par stratification

L’application de la modélisation avec Cubist pour les différentes strates géographiques suit l’algorithme suivant :

Pour chaque strate :
  1. Définition des paramètres de modélisation
    - Validation croisée sur les différents paramètres de modélisation (10 fois)
    - Sélection des meilleurs paramètres
  2. Application du modèle avec les meilleurs paramètres (100 fois)
    - Sélection aléatoire de 80% du jeu de données pour l'apprentissage et 20% pour la validation
    - Construction du modèle sur 80% du jeu de données
    - Calcul de l'importance des variables dans le modèle
    - Prédiction sur 20% du jeu de données et comparaison avec le jeu de données initial
    - Calcul des indicateurs de qualité (R2, MSE, RMSE)
  - Calcul de la moyenne des indicateurs de qualité sur le nombre de répétition
  - Calcul de l''importance moyenne des variables pour le nombre de répétition
Fin

5.2.1 Stratification par type de climat

Sur les 3 strates climatiques modélisées, les tests de modélisation présentent des r2 de

Le rôle des variables explicatives est différent. Pour les zones 3 et 4, le taux d’argile joue un rôle prépondérant et les variables

Dans la zone climatique 5, les variables climatiques sont nettement plus contributrices du modèle que les variables d’occupation du sol.

Moyenne des indicateurs de la validation croisée
r2	MSE	RMSE
0.4331	12.4	3.436
0.5129	7.843	2.79
0.687	3.921	1.972
0.6148	6.877	2.599
0.8351	7.591	2.744

TableGrob (2 x 1) “arrange”: 2 grobs z cells name grob 1 1 (1-1,1-1) arrange gtable[arrange] 2 2 (2-2,1-1) arrange gtable[guide-box]

5.2.2 Stratification par grande région d’élevage

Ci-dessous, commenter les résultats de la stratification avec les principales régions d’élevage

Moyenne des indicateurs de la validation croisée
	r2	MSE	RMSE
A	0.5476	5.272	2.269
B1	0.8325	4.482	2.105
B2	0.6509	3.633	1.889
C1	0.8691	5.374	2.306
D	0.5014	11.18	3.298
E1	0.5918	5.743	2.369
E2	0.497	8.243	2.824